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Molecular dynamics simulations are obtained and analyzed to study pairing of l-hexyl-3-methylimidazolium 
and tetrafluoroborate ions in n-pentanol, in particular by evaluating the potential-of-mean-force between 
counter ions. The present molecular model and simulation accurately predicts the dissociation constant in 
comparison to experiment, and thus the behavior and magnitudes for the ion-pair pmf at molecular distances, 
even though the dielectric constant of the simulated solvent differs from the experimental value by about 
30%. A naive dielectric model does not capture molecule structural effects such as multiple conformations 
and binding geometries of the Hmim + and BF 4 ~ ion-pairs. Mobilities identify multiple time-scale effects in 
the autocorrelation of the random forces on the ions, and specifically a slow, exponential time-decay of those 
long-ranged forces associated here with dielectric friction effects. 



I. INTRODUCTION 

Ion-pair encounter, binding, and dissociation is cen- 
tral to solution chemistry, including chemistry in organic 
solvents^ and to understanding specific electrolyte solu- 
tions in a wide range of practical settingsP We recently 
obtained an experimental determination of ion-pairing 
specifically of l-hexyl-3-methylimidazolium tetrafluorob- 
orate in n-pentanolP Those results encouraged us to pur- 
sue detailed testing of the molecular-scale description of 
ion pairing and dynamics for that system. We report 
results of that testing here. We report, among other 
results, evaluation of memory functions for ion mobil- 
ity that provide a direct signature of long-ranged ion- 
solution interactions^ a signature that we had not antic- 
ipated in advance of these simulations. 

The theoretical testing naturally relies heavily on 
molecular simulations that have become accessible in re- 
cent years. [Details of the present simulation calcula- 
tions are provided in the Appendix.] The targets of 
these calculations were the ion-pair potential-of-mean- 
force (pmf) , the ion-pairing dissociation constant implied 
by that pmf, then further the mean-square-displacement 
of the ions individually, their velocity autocorrelation 
functions, and the corresponding memory functions. To 
assist the interpretation of these quantities, we also eval- 
uated the dielectric constant of the solvent in the present 
simulation model. 



II. RESULTS AND DISCUSSION 

The analysis requires that we choose a center for the 
molecular ions of interest. Radial distribution functions 
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associated with the chosen centers, (Fig. [I]) mid-C for the 
l-hexyl-3-methylimidazolium (Hmim + ) ion and B for the 
tetrafluoroborate (BF 4 _ ) ion, characterize ion solvation. 

The CB-pair pmf, w(r), (Fig. [2]) shows strong asso- 
ciation. Comparison with a naive dielectric continuum 
model, w(r) ~ —q 2 /4irer, highlights the molecular struc- 
ture of the simulation result. The local minimum for 
r ss 0.65 nm reflects the binding of the BF 4 ~ ion to the 
ring center of the Hmim + ion. Thus multiple binding 
geometries can be a specific complication in pairing of 
molecular ions. Additionally, the ionic electric charge is 
significantly distributed over these molecular structures, 
so this dielectric model (Fig. [2]) is indeed naive for that 
reason also. 

At the longest range here, the computed pmf is fore- 
shortened due to the periodic boundary conditions, i.e., 
the mean forces between the ions along a cartesian axis 
will be zero at a boundary face for the simulated result 
but not for the dielectric model, which is expected to be 
correct at the longest range. Recognizing that distinc- 
tion, the dielectric model utilizing the dielectric constant 
evaluated for the simulation model n-pentanol solvent 
(Fig. [3| over-estimates the maximum binding free en- 
ergy of this system. This naive dielectric model predicts 
a more accurate maximum binding free energy when the 
experimental dielectric constant is utilized, but that com- 
parison has no physical significance. 

The dissociation equilibrium ratio is defined by refer- 
ence to the equilibrium 

Hmim • BF 4 *=; Hmim+ + BF 4 " , (1) 

and then 

Ki _ (PHmim+) (PBF 4 ~) /gs 
PHmimBF 4 

with p a the number density of species a. The species in- 
dicated in Eq. ([2| are identified by determining whether 
the mid-C . . . B atom pairs are within a radius r (paired) 
or not (unpaired). Basic statistical thermodynamics 
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FIG. 1. Radial distribution functions of n-pentanol O-atom 
from the ion centers (upper) middle C-atom of the imida- 
zolium ring and (lower) the B-atom of BF 4 ~, together with 
the distance-ordered decompositions (grey). For the upper 
panel n = 1, 2, 3 fills out the principal peak and n=\ partic- 
ipates in both first and second shell. For the lower panel n = 
1, . . . , 5 fills out the principal peak and n—6 participates in 
both first and second shell. 



identifies this raticP^ as 



1_ p(n = 1) 
K d p{n = 0)p B F 4 



(3) 



and we acknowledge that /?bf 4 - = /°Hmim+- P( n ) is the 
probability of observing n mid-C atoms of Hmim + ions 
within r of the B-atom of a distinguished BF4~ ion. For 
the infinite dilution circumstances of the present study, 
the ratio of probabilities (Eq. |3|) is given precisely by 
the Poisson distribution formul™ 
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4tt / e^q>[-wir')/k B T}r' 2 dr' .(4) 
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FIG. 2. The CB (Fig. [TJ pmf in n-pentanol as a function 
of the CB radial displacement. For the modelled n-pentanol 
solvent, the static dielectric constants are e/eo = 7.1 and 10.0 
at p = 1 atm and T = 298.15 K and 348.15 K, respectively. 
The dotted curve utilizes the experimental dielectric constant 
(15)P 



If the lengths on the right of Eq. ^ are in nm, and 
if Kd is mol/dm 3 , the factor for conversion of units is 
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FIG. 3. Computed dielectric constant for the simulated n- 
pentanol liquid at T — 298. 15K and p—1 atm. The shaded 
band indicates the 95% confidence interval [9.6,10.3], estab- 
lished by bootstrap resampling. See Sec. |VB| 



6.023xl0 23 /10 24 = 0.60 23 dm 3 /(mol nm 3 ), so 

1 



Ka 



(0.6023) x 4tt /J exp [-w(r')/k B T] r' 2 dr' 



(5) 



The predicted dissociation constant (Fig. |4| opera- 
tionally plateaus for r > 0.65 nm, and closely agrees with 
the experimental result. The ripple near r ~ 0.6 nm re- 
flects the second (outer or ring binding) geometry iden- 
tified with Fig. [2j and thus that second binding mode 
affects the predicted Overall, we conclude that the 
present molecular model and simulation predicts correct 
behavior and magnitudes for the pmf at molecular dis- 
tances, even though that pmf is foreshortened at long- 
range as noted above. 

The kinetics associated with the mobilities of these 
ions were characterized first on the basis of the observed 
mean-square-displacements (Fig. [5]) 

' dt ' =2 jf (v(0)-v(r))dr , (6) 

with the underlying velocity autocorrelation functions 
(Fig.© 

C(t) = (v(0).v(t))/(\v\ 2 ) . (7) 

The indicated velocities are those of the center-of-mass 
of the extended molecular ions. We also considered a 
memory function j(t) defined bjDSl 

rt 



dC(t) 



7(t - t)C{t)&t 



(8) 



with m the mass of the ion. j(t) provides the auto- 
correlation of the random forces on the ions. We em- 
phasize the connection to the forces with the notation 
Q 2 = (F 2 ) /3mk B T so that 

7 (0) = mfl 2 . (9) 
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FIG. 4. Evaluation of the dissociation constant of Eq. Q 
and comparison with experiment™ The defining integral oper- 
ationally plateaus for r > 0.65 nm, and closely agrees with the 
experimental result. The ripple near r « 0.6 nm reflects the 
second (outer or ring binding) geometry identified with Fig. [2] 
and thus that second binding mode affects the predicted 



We extracted 7(f) from C(t) on the basis of the relation^ 
e~ s Ht)dt = 7(5) = m (J- - s^j , (fO) 

which follows from Eq. with C(s) the Laplace trans- 
form of C(t). 

We utilized the Stehfest algorithm 11 to invert the 
Laplace transform numerically. 7(f) is found (Fig. uh 
to be non-negative in these cases. Furthermore, 7(F) 
displays two different time regimes, and specifically a 
slow, exponential time-decay for the longest times an- 
alyzed here. The long-time decay recalls the discussion 
of Wolynes many years agc^ of the effects of dielectric 
friction on ion mobilities. This behavior seems not to 
be have been observed until now, but Annapureddy and 
Dang have obtained the analogous result for the auto- 
correlation of the forces on stationary alkali metal ions 
in waterP^l 

III. CONCLUSION 

The present molecular model and simulation accu- 
rately predicts the dissociation constant of 1-hexyl- 
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L>bf 4 - = 2.0 x 10" 6 cm 2 /s 

^Hmim + = 1.5 X lO" 6 Cm 2 /S 
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FIG. 5. Mean-square-displacements of the center-of-mass 
of these molecular ions observed in simulations of the in- 
dividual ions in n-pentanol. The fainter peaked curves in 
the background are the time-derivatives of the mean-square- 
displacements obtained from Eq. Q. The experimental 
valu<P for the average diffusivity is 1.87xl0" 6 cm 2 /s. 



3-methylimidazolium and tetrafluoroborate ions in n- 
pentanol in comparison to experiment, and thus the be- 
havior and magnitudes for the ion-pair pmf at molecu- 
lar distances, even though the dielectric constant of the 
simulated solvent differs from the experimental value by 
about 30%. A naive dielectric model does not capture 
molecule structural effects such as multiple conforma- 
tions and binding geometries of the Hmim + and BF4~ 
ion-pairs. Mobilities identify multiple time-scale effects 
in the autocorrelation of the random forces on the ions, 
and specifically the slow, exponential time-decay of those 
long-ranged forces associated here with dielectric friction 
effects. 
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V. APPENDIX: METHODS 

Potential parameters of ionic liquid (IL) l-hexyl-3- 
methylimidazolium tetrafluoroborate were taken from de 
Andradep^ which is the AMBER force field with slight 
modification. n-Pentanol force field parameters were 
taken from AMBER. 14 Partial atom charges of Hmim + 
were derived following de Andrade's procedure, while 
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FIG. 6. Autocorrelation functions of the center-of-mass ve- 
locities of these molecular ions. The time interval shown here 
is not sufficient to obtain the self-diffusion coefficients of Fig. 
[5] On this time interval these functions are qualitatively sim- 
ilar even though the molecular structures of these ions are 
qualitatively different. 



those of BF4 - were taken directly from their workP^I 
Partial atom charges used in this work for n-pentanol 
were those developed by Kuhn.^ 

Our systems were simulated using the AMBER 10 
package in an isothermal-isobaric ensemble (NPT) with 
periodic boundary conditions. The cutoff for nonbonded 
interactions was 1.7 nm, and Ewald summation was used 
to calculate electrostatic interactions. The temperature 
was regulated with Langevin dynamics, while pressure 
was controlled by Berendsen's weak coupling algorithm. 
All C-H bonds were constrained by SHAKE algorithm.^ 
For the simulation of IL pair in n-pentanol, an ion pair 
was first equilibrated without cutoff in vacuum for 10 ns 
to get an optimized geometry Then it was placed in the 
center of a (4.8 nm) 3 cubic box, which was packed uni- 
formly with 620 n-pentanol molecules using Packmol.^ 
For the single ion/n-pentanol systems, Hmim + or BF4T 
was simply wrapped up with 620 n-pentanol molecules in 
a cubic box of the same size. Aging was 2 ns at 298.15 K 
under 1 atm with 1 fs integration time step. Then 8 ns 
production equilibrium run was performed. For the single 
ion/n-pcntanol systems, the production runs extended to 
10 ns. Configurations were saved every 1 ps for further 
analysis. For the velocity autocorrelation function and 
the mean square displacement of each ion, 1 ns trajecto- 
ries, with a time step of 2 fs, were obtained, saving the 
phase point at each 10 fs. 
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B. ra-Pentanol Dielectric Constant 

Assessment of the static dielectric constant of liquid 
n-pentanol is required in considering the molecular pair 
potential of the mean forces for Hmim + . . . BF4 - in n- 
pentanol. The dielectric constant for th is mo del was ob- 
tained by standard simulation methods! 19 * 20 * The calcu- 
lation treated 620 n-pentanol molecules under standard 
periodic boundary conditions. The results (e/eo = 7.1 
and 10.0) were extracted from averaging over 20 ns of 
simulation trajectory under constant pressure conditions, 
and at constant temperatures of 298.15 K and 348.15 K, 
respectively. These computed values for the dielectric 
constant of the modelled solvent are about 30% smaller 
than experimental values (Fig. [3]) . 
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FIG. 7. Logarithm of the normalized autocorrelation function 
of the random forces on these ions, obtained with the Ste- 
hfest algorithm^ from the velocity autocorrelation functions 
(Eq. (pi). That 7(f) is non-negative here is obvious but spe- 
cial. j[t) decays approximately exponentially for the largest 
times shown. This long-time behavior was suggested many 
years ago in considering dielectric friction on ion mobilities.™ 
This behavior seems not to be have been observed until now, 
but Annapureddy and Dang have also obtained the analogous 
result for the autocorrelation of the forces on stationary alkali 
metal ions in water.^ 



A. The pmf and WHAM calculations 



To obtain the ion-pair pmf, a reaction coordinate 
was defined as the distance between the central car- 
bon of the imidazolium ring (Fig. [IJ and boron atom 
of the BF 4 ~ ion. Window calculations then utilized a 
harmonic stratifying potential with a force constant of 
100 kcal/(mol-nm) for windows from 0.3 nm to 2.35 nm 
(at 298.15K), and from 0.3 nm to 1.95 nm (at 348. 15K), 
with an increment of 0.05 nm. Typically, the system was 
equilibrated for 1 ns to initiate each window simulation, 
and then followed by 10 ns production run. For window 
separations greater than 1.7 nm, the initial configuration 
for the next window was taken from the last configura- 
tion from the previous MD simulation. The weighted his- 
togram analysis method (WHAM]P^ was used to synthe- 
size the final pmf profile. For the calculation at 348.15 K, 
the system consisted of an IL pair and 340 n-pentanol 
molecules in a (4 nm) 3 cubic cell. 
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